# 灰色模型预测

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# 读取数据
data = pd.read_excel('data.xlsx')
data = np.array(data['A'])
lens = len(data)  # 获取数据长度 === 数据量

# 数据检验
## 计算级比
cul = []  # 存储计算级比
for i in range(1, lens):  # 计算lens-1的计算级比
    cul.append(data[i-1]/data[i])  # 计算级比 = 前一个数据的值/当前数据的值
## 计算区间
X_min = np.e**(-2/(lens+1))  # 自然底数e
X_max = np.e**(2/(lens+1))
## 检验
### 先假设检验为True
ok = True
# 对每一个计算级比的值进行范围确定
for item in cul:
    # 如果遇到一个超出范围的就设置检验为Flase
    if (item < X_min or item > X_max):
        ok = False
# 处理结束后输出结果
if (ok == False):
    print('该数据未通过检验')
else:
    print('该数据通过检验')

# 构建灰色模型GM(1,1)
## 累加数列
data_sum = []  # 储存累加数列
data_sum.append(data[0])  # 先将初始数据放入列表
for i in range(1, lens):  # 做lens-1次循环
    data_sum.append(data_sum[i-1]+data[i])  # 计算累加值
## 灰导数及临值生成数列
ds = []
zs = []
for i in range(1, lens):
    ds.append(data[i])
    zs.append(-1/2*(data_sum[i-1]+data_sum[i]))  # 根据灰色倒数公式放入数组
## 求a、b_msg
B = np.array(zs).reshape(lens-1, 1)  # 生成B矩阵
one = np.ones(lens-1)  # 生成单位矩阵
B = np.c_[B, one]  # 加上一列1 np.c_为按行接入矩阵
Y = np.array(ds).reshape(lens-1, 1)
a_msg, b_msg = np.dot(np.dot(np.linalg.inv(np.dot(B.T, B)), B.T), Y) # 根据公式对矩阵操作转置等
print('a_msg='+str(a_msg))  # 打印结果
print('b_msg='+str(b_msg))

# 对未来值进行预测
def func(k):
    c = b_msg/a_msg
    return (data[0]-c)*(np.e**(-a_msg*k))+c
data_sum_hat = []  # 累加预测值
data_origin_hat = []  # 原始预测值
data_sum_hat.append(func(0))
data_origin_hat.append(data_sum_hat[0])
'''
在这里设置向后预测，前面作为训练计算
'''
for i in range(1, lens+15):  # 多预测5次
    data_sum_hat.append(func(i))
    data_origin_hat.append(data_sum_hat[i]-data_sum_hat[i-1])
print('预测值为：')
for i in data_origin_hat:
    print(i)

# 模型检验
## 预测结果方差
data_h = np.array(data_origin_hat[0:7]).T
sum_h = data_h.sum()
mean_h = sum_h/lens
S1 = np.sum((data_h-mean_h)**2)/lens
## 残差方差
e = data[0:7] - data_h
sum_e = e.sum()
mean_e = sum_e/lens
S2 = np.sum((e-mean_e)**2)/lens
## 后验差比
C = S2/S1
## 结果
if (C <= 0.35):
    print('1级，效果好')
elif (C <= 0.5 and C >= 0.35):
    print('2级，效果合格')
elif (C <= 0.65 and C >= 0.5):
    print('3级，效果勉强')
else:
    print('4级，效果不合格')

# 画图
plt.figure(figsize=(9, 4), dpi=100)
# 这里要根据预测数量及实际数量进行修改x
x1 = np.linspace(1, 20, 20)
x2 = np.linspace(1, 35, 35)
plt.subplot(121)
plt.title('x^0')
plt.plot(x2, data_origin_hat, 'r--', marker='*')
plt.scatter(x1, data, marker='^')
plt.subplot(122)
plt.title('x^1')
plt.plot(x2, data_sum_hat, 'r--', marker='*')
plt.scatter(x1, data_sum, marker='^')
plt.show()
